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Control systems theory concerns the study of open dynamical systems that process time-dependent input 
signals (stimuli, ligands, controls, forcing functions, test signals) into output signals (responses, mea- 
surements, read-outs, reporters). Such input/output (i/o) systems may be studied by themselves, or as 
components (subsystems, modules) of larger systems. Two extreme approaches to modeling systems (or 
w , particular subsystems) are as follows Q|2]|: 

(-H ] • The black-box, purely input/output, or phenomenological description ignores all mechanistic detail 

and characterizes behavior solely in terms of behavioral or stimulus-response data. This data may be 
summarized for example by transfer functions (for linear systems and certain types of "bilinear" sys- 
tems), Volterra, generating series, or Wiener expansions (for general nonlinear systems), statistical 
models that correlate inputs and outputs, neural network and other approximate models, or machine 
learning prediction formalisms. 

\C> • • The state-space, internal, or mechanistic description, in which all relevant players ("state variables" 

^ \ such as proteins, mRNA, and metabolites) are specified, and the forms for all reactions and reaction 

Q ' constants are described. 

\D 

O ' In between these two extremes, one finds combinations ("gray box" view). Available mechanistic in- 

^ . formation is incorporated into the model, and this is supplemented by i/o data, which serves to impose 

\ constraints on (or, in a Bayesian approach, to provide prior distributions of) unknown paiameters to be 

further identified. In engineering applications, black-box modeling is common, for example, in chemical 
process control, and state-space modeling is often used in mechanical and aerospace engineering. Many 
biological models involve combinations of both approaches. 
The subfield of realization theory ("reverse-engineering") deals with the question of how much of the in- 
■ ternal system can be deduced from its input/output or black-box behavior One may pose questions about 

the structure itself as well as about parameters, once that a concrete systems structure has been hypoth- 
esized. The basic results classify symmetry groups of possible internal realizations of a given behavior, 
typically through linear algebra, algebraic-geometric, and differential-geometric techniques, depending on 
the class of systems being studied, or provide conditions for paiameter "identifiability" from i/o data. (The 
term is not to be confused with identification theory, the field that incorporates stochastic aspects so as to 
deal with noisy data. The underlying theoretical questions of what is ultimately achievable, even in ideal 
noise-free conditions, are easier to understand in a deterministic setting.) 

Key concepts used in realization theory are observability (how much information about internal dy- 
namic variables, such as chemical concentrations, can be ultimately obtained from time-varying stimu- 
lus/response measurements?) and reachability (or controllability) (what state configurations can be at- 
tained, from a given initial state by manipulating inputs, also called "open-loop" controls?). Minimal 
models explaining observed input/output behavior typically are characterized in terms of reachability and 
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observability. Also, identifiability of parameters can be reduced to an observability problem. Many tests 
exist for observability and controllability, for linear and more generally for large classes of nonlinear sys- 
tems yj, although the theory of nonlinear controllability is still very incomplete. 

There are many other important subfields of control theory, among them the study of disturbances and 
robustness (how does one quantify the effect of actuator and sensor disturbances, and the robustness to 
unstructured perturbations?), feedback (how does one design feedback laws or "closed-loop controls" so 
as to enhance desirable system characteristics such as stability and robustness to structural perturbations, 
noise, and other disturbances?), optimal control (among those open-loops controls which serve to drive one 
internal configuration of state variables to another, which ones are optimal in appropriate senses?), and the 
question of hybrid or discrete/continuous interfaces (how do digital sensors, actuators, limited-bandwidth 
communication channels, and feedback devices interact with systems described by analog real-valued 
signals?). Each of these subfields brings up its own set of challenges special to systems biology [3^ 4J- 

Here, we make some elementary remarks regarding how the choice of inputs being applied to a system 
(pulses, steps, etc) affects parameter identifiability. To keep the exposition accessible, we work out some 
extremely simple explicit examples, instead of discussing general theoretical matters. Some of the mate- 
rial, such as that covering linear systems, is standai^d |1|. The nonlinear examples illustrate ideas from a 
paper in preparation. 

Various Input Classes 

We discuss next how the choice of input signals impacts the theoretical possibility of identifying system 
responses (black-box behavior) and/or internal structure (state-space view). 
Normally, the analysis of responses to small signals (or small perturbations from steady states, leading to 
"weakly activated" systems) can be handled by linear control theory, which provides a rich toolkit for the 
analysis of biochemical networks. On the other hand, if relevant signals are large, or if nonlinear effects, 
such as receptor desensitization or dimerization, affect signaling at a comparable time-scale, then linear 
tools do not suffice and new techniques must be brought to bear. Nonlinear control systems analysis has 
undergone an exceptional period of development and maturation within the last 20 years, largely driven 
by application needs from aerospace, mechanical, and other engineering disciplines. The tools that have 
been developed reflect the requirements from these fields, hence the emphasis on studying certain classes 
of systems, such as Hamiltonian dynamics, which are less relevant to biochemical networks. Although 
many of the fundamental ideas are still applicable to biology, much further theoretical work is needed. 

Small-signal analysis: equivalence of step, impulse, and pulse responses 

When probing systems with signals of small magnitude, the full power of the well-developed linear control 
theory can be taken advantage of. To illustrate the problems that arise in this context, we use a very simple 
example, that of a single species X whose concentration x{t) is described by an ordinary differential 
equation: 

^{t) = -ax{t) + hu{t) 

in which the variable u{t) measures the, generally time varying, concentration of an external input U, and 
a, b are unknown positive constants. (From now on, we will drop the time arguments t when clear from 
the context.) 

For example, X might represent the concentration of bound receptors, and U the concentration of a 
ligand which is present only in very small quantities. In this case, a is the degradation rate of X, due 
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to internalization, ligand unbinding, and other effects, and b represents the concentration of unoccupied 
receptors, which we consider as constant by assuming that their abundance is large compared to the ligand. 
(If the concentration is not constant, then a linear model is not appropriate.) Generally, variables such as x 
and u represent not absolute concentrations but rather perturbations from a steady state. 
We assume that the reporter being used provides a measurement y{t) = cx{t) equal to the product of the 
concentration of X and some (unknown) positive constant c. Mechanistically, the form of y might arise 
from an additional differential equation: 

dy 

£— = —y + cx 
dt ^ 

which models the formation and decay of the reporter Y (we took a unit decay to simplify) and < e <C 1 
indicates that this process happens at a much faster time-scale than the production of x. Using a quasi- 
steady state approximation (singular perturbation), we set e = and consider y = cx as the output. 
Another way in which an output of this type could arise is as a model of unknown measurement device 
calibration. 

For any input u{t) for t > 0, and x(0) = 0, we have, for t > 0: 

y(t) = {k * u)(t) = / k{t — r)u{r)dr 
Jo 

(star denotes convolution), where 

k{r) = ce""*6 . 

The function k{t) is called the impulse-response of the system, because a unit impulse u = Sq at time 
(or, in practice, a very narrow pulse with unit area) elicits the response y{t) = k{t). The function K{t) is 
called the step-response, because for any step input u{t) = uq for t > 0, the output is y{t) = K{t)uo for 
t > 0. 

Knowledge of K is equivalent to knowledge ofk, because k = K' and K = Jq k. Therefore, knowing the 
time-series response to any single step input is sufficient to determine K{t) (just divide: K{t) = y{t)/uo), 
which in turn (convolution formula) determines the response to any other possible input, for instance 
pulses. 

In fact, knowing the response to basically any single input is theoretically enough for finding k. A proof 
using Laplace transforms, is clear from the fact that y{s) = W{s)u{s), where W{s) is the frequency 
response of the system, that is to say, the Laplace transform of k. Thus, as long as u{s) does not vanish, 
one can solve for = y/u, and k can then be obtained by Laplace inversion. In practice, of course, 
one does not perform just a single input/output experiment, since measurements will be noisy, nor does 
one obtain k as K', since differentiation of an estimated step function is numerically difficult. Several 
measurements are averaged, through a least-squares regression or another technique. This is discussed at 
length in textbooks in identification theory. It is easier, however, to explain the theory in the noise-free 
case. 

Steady state measurements do not provide enough information about transient behavior 

Steady-state experiments measure the steady state output y^s of the reporter variable as a function of the 
magnitude of a step input. With the above example, yss = cxss = ^^o- In other words, the only infor- 
mation about the three parameters a, 6, c that can be obtained from such experiments is the combination 

7 = he/ a (called the steady state gain of the system). For example, the two systems (i) dx/dt = —x + u 
with y = X and (ii) dx/dt = —2x + u with y = 2x have the same steady-state gain, 7 = 1. The i/o behav- 
iors of these two systems are different, however: when using a step input u{t) = 1, and initial condition 
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x{0) = 0, the output of the first system is y{t) = I — e *, but the output of the second system is a different 
function, y{t) = 1 — e^^*. Only transient measurements reveal this distinction. 

Parameter unidentifiability leads to equivalence classes of parameters 

Still with the same example, we note that the product be is an i/o invariant of the system, because b and 
c appear together as a product in the convolution formula. No matter what input is applied (assuming 
x{0) = 0), two systems for which the same product be is the same will give the same output. No amount 
of i/o experimentation will permit finding b and c individually. 
On the other hand, the decay a can be determined uniquely from i/o experiments. For example, with a unit 
step input, we may measure y{t) = K{t), from which we may evaluate k{t); hence k{0) = K'{0) = eb 
and A;'(0) = K"{0) = —eab are known, and a = —K"{0)/K'{0) can be estimated. One says that the 
two parameter triples (ai, bi,ei) and (02, &2, C2) are equivalent if they cannot be distinguished by any i/o 
experiment. Thus, we showed that two triples are equivalent if and only if biei = 62C2 and ai = 02- 
Another way to say state this fact is through a one-parameter group of symmetries: (a, b, c) is equivalent 
to (a, Tb, eT^^) for all T / 0. (These symmetries amount to rescaling the concentration x{t): in different 
units, one has a new variable z{t) = Tx{t), which satisfies the differential equation dz/dt = —az + Tbu, 
and the reporter is now and y = eT~^x.) 

A far-reaching generalization of this result is available for linear systems of arbitrary dimension n (num- 
ber of species), not merely n = 1 as in our example: systems dx/dt = Ax + bu,y = ex, where now A is 
an n X n matrix, 6 is a column n-vector, and c is row n-vector. Under genericity conditions ("reachability 
and observability"), the following theorem is true 1 1 1: two triples {Ai, 6j, q), i = 1, 2 are i/o equivalent if 
and only if {TAiT-'^,Tbi,cT{'^) = {A2, 62, C2) for some invertible matrix T. (The theorem is also true 
for multivariable inputs and outputs, in which case both B and C are matrices.) 

Continuing with our example, we may think of our system as a "gray box" in the sense that we do not 
know the complete mechanism (all the parameters) but neither do we have complete ignorance, since we 
are able to measure a and the product be. Now suppose that a new experimental design allows us not merely 
to measure y{t) but also do directly measure the pai^ameter c (for example, the binding of X to species Y 
may be characterized by means of a biochemical experiment). Now it is possible to also determine b, since 
b = K'{Q) je. Or perhaps, instead of a biochemical experiment, a new probe is devised, which allows one 
to directly measure the concentration x{t). Then, knowing that dx/dt = —ax + bu, and from knowledge 
of x{t) and x'{t) at two values t = ti, one can determine both a and b, and hence, c = K'{0)/b as well. 
In either case, we have gone from a "gray box" description to a complete mechanistic description. Partial 
knowledge of 6 or c can also be understood in these terms: suppose that we know that K'{0) = 1 and that 
an imprecise biochemical assay provides the range c € [0.01,0.1]. Then, we know that b G [10, 100], a 
"less gray" level of description than total ignorance of b. 

These ideas are particularly useful in a Bayesian context. The i/o knowledge given by the value of a 
and the product be, combined with physically realistic ranges for the unknown b and c, can be used to 
provide priors for parameter estimation through further experiments (or to evaluate the reliability of using 
parameters from a homologous protein, for example). That is, the invariants provide a "prior belief" about 
the parameters vr = (a, b, c), summarized as a probability -P(vr). After an experiment e, this probability 
may then be updated to obtain a posterior distribution -P(vr|e) through an apphcation of Bayes Rule: 

^ ' ^ P{e) ' 

where -P(e|7r) is the likelihood of the experimental results under the assumed parameters tt. 
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Step inputs are not sufficient for nonlinear identification 



There are strong theoretical justifications for the need for rich input classes in the context of reverse en- 
gineering biological systems. We use a simplified, and hence somewhat unrealistic, example in order to 
explicitly illustrate this point, but the phenomenon scales to larger systems. 
Consider the following system, which describes the time-evolution of the concentrations x{t) and z{t) of 
two chemical species X and Z: 



dx 

Itt 
dz 

'dt 



-Xx + 
-Xz + u 



with initial concentrations ,x(0) = and z(0) = 0, and where ti{t) is the (possibly time-varying) con- 
centration of an input (for instance a ligand) U that helps promote formation of X and Z. Our goal is to 
estimate the unknown parameter A, the degradation rate of X and Z. Note that as long as no input U is 
appUed {u{t) = 0), the concentrations of X and Z remain at steady state (zero). 

We assume that the only available measurement is the (time- varying) concentration y{t) of a reporter Y, 
which we take as a function of the instantaneous concentrations of X and Z according to the formula: 

y = X + u{a — z) 

were a is some known positive constant. (One could view y as obtained through a quasi-steady state 
approximation applied to 

dy , , 

£— = —y-\-x + u{a — z) 
dt 

with < e 1, which describes a reaction in which a certain substance A dimerizes with U in order to 
promote formation of y, but A is sequestered by Z with 1-1 stoichiometry, so the free amount of A equals 
a — z = ^totai ~ Abound, bccausc the concentration of Abound equals the concentration of Z.) 

If we use a step function input, u{t) = uq for t > 0, then x{t) = z{t)uQ, and therefore y{t) = auQ, 
independently of the actual value of A. Thus, with step inputs, no matter how many experiments (different 
values of uq) are carried out, it is impossible to obtain any information about X. 

On the other hand, suppose that we use a unit ramp: u{t) = t. Then y""{t) = —X{Xt — 2)e~^^, and so 

\im y""(t) = 2 A 

and therefore A can be identified by differentiation. (Of course, differentiation will not be used in practice, 

nor is it needed since one may use other techniques for estimating A, but we are just proving theoretical 
identifiability.) In summary, a single ramp input is enough to identify A, but no amount of even complete 
transient time-measurements of step responses will achieve this result. 

Instead of ramps, we could also use pulses. Suppose that we consider a unit pulse: u{t) = 1 for < t < 1 
and u{t) = for t > 1. Now, x{t) = z{t) = \{l - e"^*) for t G [0, 1] and x{t) = z{t) = \{e^ - l)e-^* 
for t > 1, so y{t) = x{t) + 0{A — z{t)) = x{t) for t > 1, which means that 



lim y'{t) = e 



-X 



and therefore A can be calculated as: 



A = - In ^1 + limy'(t)^ . 
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Once again, a non-step signal is successful in identifying the parameter A, but steps are not enough. 
Similar examples can be given with more parameters and more complicated systems. Also, we took the 
same degradation rates for X and Z, because we are considering a noise-free ideal situation. If these 
degradation rates would not necessarily be equal but are close, then y{t) will depend on A, even in the 
step case. However, the dependence will be through a very small multiplier coefficient, indicating a very 
small sensitivity and therefore very high estimator variance (Cramer-Rao bounds), thus rendering estima- 
tion practically impossible. It is possible to prove (paper in preparation) that for bilinear systems (those 
in which the only possible products are between inputs and states), pulses provide as much theoretical 
information as any possible time-varying input signals. 

Remarks 

We have only discussed parameter identification, as opposed to structure identification. For small-signal 
(linear) analysis, the two problems are the same (at least if an upper bound on the number of species in 
the system is known), since zero entries in the differential equations indicate non-existent interactions. 
For nonlinear systems, the problem of identification of structure is far more difficult. It is the subject (at 
least for the ideal noise-free case) of nonlinear realization theory |2|. The parameter estimation prob- 
lem is a special case of a nonlinear observability problem (thinking of unknown parameters as constant 
states), and a rich theory allows the use of Lie-algebra and differential-algebra techniques for determining 
observabiUty d- 

There are results available characterizing the number of experiments required for identifying a given 
number of parameters in a system |5|. Another relevant theoretical mathematical fact is that generic 
"random" single inputs are enough for identification t2J, and this fact is very useful in the context of 
identification of certain biochemical systems l6l . 
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